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IMAGING WITH POWER CONTROLLED SOURCE PAIRS 


PATRICK BARDSLEY AND FERNANDO GUEVARA VASQUEZ 


Abstract. Scatterers in a homogeneous medium are imaged by probing the 
medium with two point sources of waves modulated by correlated signals and 
by measuring only intensities at one single receiver. For appropriately chosen 
source pairs, we show that full waveform array measurements can be recovered 
from such intensity measurements by solving a linear least squares problem. 
The least squares solution can be used to image with Kirchhoff migration, even 
if the solution is determined only up to a known one-dimensional nullspace. 
The same imaging strategy can be used when the medium is probed with 
point sources driven by correlated Gaussian processes and autocorrelations 
are measured at a single location. Since autocorrelations are robust to noise, 
this can be used for imaging when the probing wave is drowned in background 
noise. 

Keywords. Intensity-only imaging, travel-time migration, noise sources, au¬ 
tocorrelation. 

AMS Subject Classifications. 78A45, 78A46, 35R30 


1. Introduction 

Scatterers in a homogeneous medium can be imaged by probing the medium with 
a wave emanating from a point source, and recording the reflected waves at one 
or more receivers. An image of the scatterers can be generated by repeating this 
experiment while varying the position of the source and/or receiver and using classic 
methods such as the Kirchhoff (travel time) migration (see e.g. [I]) or MUSIC (see 
e.g. 0). We are concerned here with the case where only intensity measurements 
can be made at the receiver; destroying phase information that migration and 
MUSIC need to image. Intensity measurements occur e.g. when the response 
time of the receiver is larger than the typical wave period or when it is more cost 
effective to measure intensities than the full waveform. This is typical in e.g. optical 
coherence tomography PH [23] and radar imaging [7] . Another situation is when 
the wave sources are stochastic and the measurements consist of correlations of the 
signal recorded at different points pziEn]. In the special case of autocorrelations 
(i.e. correlating the signal with itself), the Wiener-Khinchin theorem guarantees we 
are measuring power spectra (see e.g. |18jb another form of intensity measurements. 

The setup we analyze consists of an array of sources and one single receiver that 
can only record power spectra, i.e. the intensity of the signal at certain frequency 
samples. A crucial assumption for our method is that we can use source pairs, 
meaning we can send correlated signals from two different locations. Thus we allow 
for known delays or attenuations between the signals in a source pair. In acoustics, 
one way of achieving this would be to drive two transducers in an array with the 
same signal. With light, one could use an incoherent plane wave with wavefronts 
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parallel to a configurable mask. The mask lets light through one or two small holes, 
whose locations can be controlled. 

Our method can be used for imaging from both measurements of intensities 
(11.11 and autocorrelations (11.2). 


1.1. Intensity only measurements. One way to deal with intensity measure¬ 
ments is phase retrieval, i.e. first recovering the phases from intensity measure¬ 
ments, and using this reconstructed field to image. In diffraction tomography, inten¬ 
sity measurements at two different planes can be used to recover phases [161 US IH] ■ 
If additional information is known (e.g. the support of the scatterer), intensities at 
one single plane can be used IS US US- Total or partial knowledge of the incident 
field can also be exploited to image from intensities at one single plane |S]- 

Chai et al. |3] take a compressed sensing approach to image a few point scatterers 
exactly. With knowledge of the incident field, the location of the scatterers can be 
resolved in both range and cross-range with monochromatic measurements. The 
same ideas can even be used to deal with multiple scattering [5]- Novikov et al. 
pi] use the polarization identity 4Re(it*t;) = ||m -I- — ||tt — u|p, u,v £ , 

and linear combinations of single source experiments to recover dot products of two 
single source experiments from intensity data. MUSIC can then be used to image 
with this quadratic functional of the full waveform data. 

Here we do phase retrieval assuming knowledge of the intensity of the incident 
field. Our illumination strategy using source pairs does not require direct manipula¬ 
tion of phases or addition/subtraction of wave fields. We reduce the recovery of the 
total field to a linear system with a one-dimensional nullspace which we can write 
explicitly in terms of the incident field. There is one (very sparse) linear system per 
frequency sample to solve, and the linear system has size comparable to twice the 
number of source positions. Intuitively we are recovering a field in from 2N (or 
more) real measurements. We show that vectors in the one-dimensional nullspace 
do not affect Kirchhoff migration. Hence we can use, without modification, Kirch- 
hoff migration and its standard range and cross-range resolution estimates (see e.g. 

my 

1.2. Correlation based methods. In seismic imaging, correlations of traces (or 
recordings) at many receivers have been used to image the earth’s subsurface, es¬ 
pecially when the wave sources and their locations are not well known |25l [271126]. 
The idea is that correlations of the signals at two different locations contain infor¬ 
mation about the Green’s function between the two locations, and this information 
can be exploited to image the medium and any scatterers. This principle can even 
be exploited to do opportunistic imaging with ambient noise nnuniiii]. Cross¬ 
correlations can also be used to image scatterers in a random medium [allllllls]. In 
radar imaging, the measurements are in fact correlations |7] , and so even stochastic 
processes can be used instead of deterministic signals pa ISO]. 

The method we present here can also be used to image scatterers using autocor¬ 
relations. We show it is possible to form an image by exploiting angular diversity 
in source pairs instead of cross-correlations among different receivers. Just as in 
the intensity measurements case, we are able to recover (up to a one dimensional 
nullspace) full waveform array measurements. One advantage of using autocor¬ 
relations instead of cross-correlations is that the data acquisition at the (single) 
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receiver is simpler. The drawback is that our illumination strategy requires to il¬ 
luminate with pairs of sources, but also with each of the sources in a pair on its 
own. To get the same full waveform data as an array with N sources, we need 
at least 3fV different experiments. Another advantage of using autocorrelations is 
that the measurements are extremely robust to noise. As an example, our numer¬ 
ical experiments show that it is possible to image scatterers with an array that 
is sending noise from all possible source positions; the only assumption about the 
noise being that all the sources are independent stochastic processes except for two 
correlated sources whose positions we can control. Because of this robustness, it 
may be possible to use our imaging method in situations where the medium is to 
be probed in a non-intrusive way, i.e. active imaging with waves that are of the 
same magnitude as the ambient noise. 

1.3. Contents. The particular physical setup we consider is described in ^ The 
illumination strategy with source pairs is explained in which leads to a phase 
retrieval problem that can be formulated as a linear system (Q. The least squares 
solution to the linear system is then used as data for imaging with Kirchhoff mi¬ 
gration, and we show that this gives essentially the same images as full waveform 
data (©■ The extension to stochastic source pairs is given in ^ Then we show 
that our method is robust to additive noise when using autocorrelations (Q. Nu¬ 
merical experiments illustrating our method are given in ^ and we conclude with 
a discussion in ^ 


2. Array imaging for full waveform measurements 


Here we introduce the experimental setup we consider (^2.1) and briefly recall 
the classic Kirchhoff migration imaging method (^2.2). 


2.1. Experimental setup. The physical setup is illustrated in hgurej^ We probe 
a homogeneous medium with waves originating from N point sources with locations 
Xs G A, s = l, 2, For simplicity we consider a linear array in 2D or a square 

array in 3D, i.e. A — [—a/2, a/2]'^“^ x {0}, where c? = 2 or 3 is the dimension. Our 
imaging strategy imposes only mild restrictions on the source positions, so other 
array shapes may be considered. Waves are recorded at a single known receiver 
location x^- 

The total field generated by the array (or incident field) can be written as 

(1) Uinc{x,uj) = go{x,uj)'^f{uj), 


where 

(2) goix,u}) 


Go{x, Xi,uj),Go{x,X2,uj),..., Go{x, xn,uj) 


G C 


N 


and the source driving signals are /(w) = [/i(a;),/ 2 (w)i ■ ■ •;/Af(w)]^- Since we 
assume waves propagate through a homogeneous medium, we used the outgoing 
free space Green function, 

^ y\), 

Go{x,y,u}) = < exp[zfc|:g- y\] 

[ 47r|£-:y| 


( 3 ) 


for d = 2, 
for d = 3. 
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A 



Figure 1. Physical setup for array imaging with an array A of 
sources a?s and a single receiver Xr- The scatterer is represented 
by a compactly supported reflectivity function p{x). 


Here is the zeroth order Hankel function of the first kind, k = lo/cq is the 
wavenumber, uj is the angular frequency and Cq is a known constant background 
wave speed. For functions of time, the Fourier transform convention we use is 

(4) 7(w) = j f{t)A^*dt, and f{t) = ^ J where / S L'^{E.). 

The scatterers we want to image are represented by a compactly supported 
reflectivity function p{x). Under the weak scattering assumption (i.e. p <C 1), we 
can use the Born approximation to the total field at the receiver 

(5) u{xr,uj) = (go+p^f, 
where the array response vector is 


( 6 ) 


p{x,uj) = 


dyp{y)Go{x, y, u})go{y, w). 


2.2. Kirchhoff migration. By e.g. using illuminations f{oj) = e^, i = 1,..., 
corresponding to the canonical basis vectors, it is possible to obtain the array 
response vector p(xr,uj) from the measurements ([^. The scatterers can then be 
imaged using the Kirchhoff migration functional (see e.g. m) which for a single 
frequency uj is 

(7) rKM[p,w](y) = Go{y,Xr,uj)go{y,uj)*p{xr,uj), 

where y represents a point in the image. This image has a Rayleigh or cross-range 
(i.e. in the direction parallel to the array) resolution of AL/a, where L is the 
distance from the array to the scatterer (see e.g. m)- To get range (i.e. in the 
direction perpendicular to the array) resolution we need to integrate rKM[Pjw](y) 
for frequencies uj in some frequency band B = [—ujmax, —<-^min] U [uJmim^max], the 
same frequency band of the signals /(w). The range resolution is then ct^jiuimax — 
^min) (see e.g. d]). We discuss this imaging functional further in section [§ 


3. Intensity only measurements 


We start in §3.1| by describing a source pair illumination strategy for intensity 
measurements of the total field \u{xr,uj)\'^. With this strategy, the problem of 
recovering the array response vector p can be formulated as a linear system ((3.2). 





IMAGING WITH POWER GONTROLLED SOURCE PAIRS 


5 


3.1. Illumination strategy. The data we use comes from probing the medium 
with Np source pairs that are sending signals with known power and phase differ¬ 
ence. Since the number of distinct source pairs out of an array with N sources is 
N{N — l)/2 we must have Np < N{N — l)/2. We assume the power and phase 
differences remain the same for all Np illuminations. The case where these quan¬ 
tities depend on the source pair is left for future studies. To be more precise, the 
illumination corresponding to the m—th source pair {i{m),j{m)) G {1,... ,iV}^ is 


( 8 ) 


fm{^) — 


a{uj) 

pluj[ 


where e 


We emphasize that only \a\^, |/?P and the phase difference = arg(a/3) is 

assumed to be known for the signals a and p. A particular case is when the same 
signal is sent from the source pair, i.e. P = a and </>(w) = 0. 

The intensity of the field Um arising from the source pair illumination is 


(9) 


ixr,ui)\ = g^f^f^g = g*f„ 


aP 


Pa 




where we used g = go + P- Note that since aP = |a||/3|e*‘^, the inner 2x2 Hermit- 
ian matrix is uniquely determined by the magnitudes of a and P and their phase 
difference p. By using the single source reference illumination we additionally 
measure 


( 10 ) 


=g*e,ejg, 


ioT i = 1,... ,N. 


The data we exploit to recover p is obtained by subtracting the appropriate refer¬ 
ence illuminations (10) from ([^, that is 


— I’^ml |g^| 


I 


nsV)i 


= g 


0 aP 
Pa 0 


^ m9' 


3.2. Phase retrieval problem as a linear system. By recalling that g = go+P, 
the measurements are 


dmiXr^ ai) 


(fifo +P)*Fm 


0 

Pa 


aP 

0 


FlCfl-o +P)- 


To make the following expressions concise, we denote by D the Hermitian matrix 


( 11 ) 



aP 

0 


By the weak scattering assumption, we may neglect the quadratic terms in p and 
collect all measurements for m = 1,..., Np as a single vector d G 


di(Xr,U!) 


1 

■ 9o*FiDFT - 

\ 

d2(Xr,Uj) 

« d{xr, uj) = Re 


9 SF 2 DFT 

{go + 2p) 

_djV^(Xr,Uj)_ 


\ 

9oFApDF](f^_ 

/ 


= M(a;r,a;) 


Re (go + 2p) 

Im(go + 2p) 


( 12 ) 
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Figure 2. An example illustrating the strategy to choose the 
source pairs for = 8 source positions. Each source position 
is represented by a node in the graph, and source pairs are repre¬ 
sented by edges. The first 5 source positions are in the circle. 


where the Np x 2N real matrix 


is given by 


(13) 


w) = 


Re(gSFiDFl) 

Re(gSF2DFT) 


-Im(g5FiDF7) 

-Im(g5F2DFj) 


R-e(gSFArpDFL ) ) 


Note that by construction, the matrix M has at most 4 non-zero elements per 
row, and is thus a very sparse matrix for N large. 


4. Analysis of the phase retrieval linear system 

We now address the question of whether there is enough information in the 
measurements d G to recover the array response vector p G . The main 


result of this section is Theorem 4.3 where we show that with appropriately chosen 
pairs of sources, (i.e. the Moore Penrose pseudoinverse of M times d) gives p 
up to a complex scalar multiple of the vector Qq, which is known a priori. 

Let us first consider the case where we take measurements using all possible 
source pairs, i.e. that Np = N{N — l)/2. Clearly, we need A^ > 5 to guarantee 
that Np > 2N, i.e. that the matrix M has more rows than columns and the system 
d — M[Re(go + 2p)^,Im(go + 2p)^]^ is overdetermined. 

Instead of using all possible source pairs, we use the following strategy which for 
N > 5, guarantees Np = 2N. 

Strategy to choose source pairs: 

(1) All 10 distinct source pairs between the source positions {1,..., 5}. 

(2) For source position s > 5, choose any two different source pairs of the form 
(s,i) and (s, j) where i,j G {1,... ,5}. 

This strategy is illustrated in figure More source pairs can be added without 
affecting the recoverability of p (Theorem 4.3). We now make the following as¬ 
sumption on the first 5 source positions. 

Assumption 4.1. We assume the receiver is located at a position Xr such that for 
i,j = 1,..., 5, the vector go = go{xr, uj) satisfies 

(14) Re{go)^^0, Im{go)^^0, and Re{go) Jm{go) ■ ^ Re{go) Rm{go) ■■ 
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Additionally for one pair i, j € [ 1 ,..., 5] we assume 

cos(^) (^Re{go) .Im{go) . - Re{go) .Im{go)+ 

- siii((^) {Re{go)^Re{gQ) . + Im{go).Im{go 


(15) 


This assumption is by no means necessary for the end result (Theorem |4.3[ ) to 
hold, but it is sufficient. If d = 3, condition (14) is equivalent to the geometric 
condition 

(16) \xi — Xr\ fi- and \xi — Xr\ — \xj — Xr\ ^ for all j, j = 1,..., 5, 


while condition (15) implies for one pair i,jG [1,..., 5] that 
(17) 


\Xr — X-i \ — \Xr — Xf\ 


2 277 ^ 


Here the set {XjlAfL is the set of all integer multiples of A/2, where A = is 

the wavelength. In d = 2 , conditions similar to (16) and ([T7|) are sufficient when 


the sources and receivers are far apart because of the Hankel function asymptotic 

= \ — exp[ 7 (t — ( 7 r/ 4 ))](l + 0{l/t)), as t ^ oo. 

V 7rt 

Lemma 4.2. Provided a ^ 0, P ^ 0, ReiaP) 7 ^ 0, the source pairs are chosen with 
the a 

(18) 


the above strategy and assumption 4-1 holds, the matrix M = M(a;r,a;) satisfies 

-Im{go{Xr,uj)) 


null M = span 


Re{go(xr,uj)) 


Proof. For clarity of exposition, we adopt the notation 

Oi = Re{go)^, bi = Im(go),, 

for i = and with go = go{xr,uj). The proposed vector spanning the 

nullspace is = [—Im(go)^, R.e(go)^]^ and has components Vi = —bi and 

Wi = at for i = 1 ,... ,N . 

The proof is by induction on the number of sources N. For the purpose of 
the induction argument, we denote by the measurement matrix M(£c/,w) 

corresponding to N sources, which if we use the strategy explained above, must be 
a 2N X 27V real matrix. For the base case TV = 5 of the induction, can be 
written as 




0 

0 

0 

Bt 

Bf 

0 

0 

0 

^3 

0 

At 

0 

0 

Bt 

0 

Bf 

0 

0 


0 

0 

At 

0 

Bt 

0 

0 

Bf 

0 


0 

0 

0 


Bt 

0 

0 

0 

Bf 

0 

^ 3 - 

At 

0 

0 

0 

Bt 

Bf 

0 

0 

0 


0 

At 

0 

0 

Bt 

0 

Bf 

0 

0 

Ao 

0 

0 

At 

0 

Bt 

0 

0 

Bf 

0 

0 

A 4 


0 

0 

0 

Bt 

Bf 

0 

0 

0 

Ao 

0 


0 

0 

Bt 

0 

Bf 

0 

0 

0 


At 

0 

0 

0 

Bt 

Bf 


where we have used 

(19) = |a||/3|(cos((/)ai ± sin((/)) 6 j), Bf = |a||/5|(cos((/)6j ± sin((/)ai). 
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Using the expressions (19), the leading principal 9x9 minor of is 


''ug.ugl = -4|a|®|/3|®cos^((/)) (cos((/))( 630 i - 6103 ) + sin((/))( 63 &i 

05(6302 - 6203)(620i - a26i)(65a4 - 0564). 


0301)) X 


Therefore if assumption 4.1 holds and cosfj) 7 ^ 0 (which we get from Re(a/3) 7 ^ 0), 


we must have rankM^®) > 9. By direct calculations, we have that 

null = span I [— 61 ,..., — 65 , oi,..., 05 ]^} . 

Thus the base case TV = 5 holds and rank = 9. 

For the induction hypothesis we assume that N > 5 and that 

null = span { [-b^, , 

where a = [oi,..., and b = [ 61 ,..., 6 ^ 1 ]^- If the first 2N source pairs to 
construct are chosen in exactly the same way as the source pairs to construct 

and the last two source pairs are, e.g. (iV +1,1) and (A^+ 1, 2) we must have 
for any v,w £ and Vff+i,WN+i G K that 


( 20 ) 


Hence if uat+i,wat+i]^ G nullM(^+^\ then we must have G 

null i.e. there is some real fc 7 ^ 0 such that v = —kb and w = ka. Equating 


V 


mw 

V 


Vn +1 



w 


w 



+ B^ WN+i 



.^Af- 1 - 1^2 + A^Vn+i + + i ?2 WaT-I-I, 


the last two components of ( 20 ) to zero and using that Vi = —kbi and Wi = kai for 


i = 1 , 2 , one gets the linear system 


'At 

Bf' 


VN +1 


At 

Bf\ 


_Wn+ 1 _ 



fc^AT+l^l “ *5)v+i“i 

Since = |ap|/ 3 p(ai 62 — 7 ^ 0, the unique solution to this 

system is uat+i = —kh^+i and w^+i = fcoAr+i- Thus the desired result holds for 
any N > b. □ 

In figure]^ we show the condition number of M(a;r, w) (i.e. ai(a 2 N-i the ratio of 
the largest singular value to the smallest non-zero singular value ) over a frequency 
band. The experimental setup is that given in and corresponds to sending 
exactly the same signal from both locations in a source pair (i.e. a = f3 and (j) = 0). 
Figurej^a) shows the condition number of M with chosen so that assumption 4.1 
is satisfied, while figure [^b) shows the condition number of M with chosen so 
that assumption |4. 1 1 is violated for some frequencies. In both cases, we see improved 
conditioning by using more than 2N source pair experiments. 

We now tie to the array response vector p. 


Theorem 4.3. Under the assumptions of lemma it is possible to recover p = 
p{xr,uj) from the intensity data d up to a complex scalar multiple of Qq = go{xr,uj), 
more precisely, M^d determines the vector p + (^gQ where 

C = C{Xr,U}) = --1 - 


( 21 ) 


9590 
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(a) (b) 



Figure 3. Condition number of M , w) with receiver location x^. 


is violated for some frequencies. The number of source pair ex¬ 
periments used is Np = N{N — l)/2 (in red) and Np = 2N (in 
blue). 


chosen so that (a) assumption 4.1 is satished, (b) assumption 4.1 


d 


Proof. Recalling the form of our data we have 

Re{go + 2p) 

_Im(go + 2p) 

By lemma [T^ the matrix M has a one dimensional nullspace therefore 

-lm(go) 


MU = 


Re (go + 2p) 
Im(go + 2p) 


-c 


Re (go) 

where C € K is found by enforcing orthogonality with [—Im(go)^, Re(go)^]^, i.e. 

C = -^[-Re(go + 2p)'^lm(go) -f Im(go -h 2p)'^Re(go)] = . 

go9o go9o 

Thus from we can get the vector 

i[Re(go-h2p)-HCIm(go)]-h^[Im(go-f2p)-CRe(go)] = ^go+P-^CSo =P + Cgo, 
where the scalar ( = Ci^r, w) S C is given by (211. □ 


5. Kirchhoff migration imaging 

We now show that we can image with the reconstructed field p + (go instead of 
p by using Kirchhoff migration. This is because the Kirchhoff migration image of 
Cgo is negligible compared to the image of p for high frequencies. In order to show 
that this nullspace vector does not affect the imaging, we need to make sure the 
receiver satisfies the following condition. 

Assumption 5.1 (Geometric imaging conditions). For a scattering potential with 
support contained inside an image window W, we assume Xj. satisfies 


( 22 ) 

for s = 1,..., iV and y G W. 


X. - Xr 




Xs-y 

\x8 - y\ 


One way to guarantee assumption |5.1| h olds is to place the receiver at location 
Xr outside of the shaded region in figure 
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Figure 4. Given an array A and a region W containing the scat- 
terers to image, assumption 5.1 ensures the receiver location x^- 
is outside of the blue shaded region. This guarantees the Kirch- 
hoff images using data p and the recovered p + (go are essentially 
the same. The positive ray in the direction Xs — y for particular 
Xs & A and y G W is indicated in red. If x^ is outside the blue 
shaded region, we have (Xg — Xr)/\xs — Xr\ ^ {xg — y)/\xs — y\ 
for all Xg G A and all y G W. 


Theorem 5.2. Provided assumption \5.1\ holds, the image of the reconstructed array 
response vector is 


FkmIp + CSo, w] (y) « TxmIp, w] (y) • 

Proof. First we approximate the Kirchhoff imaging functional Q by an integral 
over the array A, i.e. 

(23) 

rKM[C9o.‘^](y) = G{Xr,y,Uj)go{y,Uj)*C{Xr,Uj)go{Xr,Uj) 

^C{Xr,uj) ( dXgC(,Xg)e^p{iu:c^Wxg-Xr\-\Xg-y\-\y-Xr\)), 

JA 


where the symbol ~ means equal up to a constant and Cixg) collects smooth 
geometric spreading terms. 

Let us first use the stationary phase method (see e.g. DP) on the integral over 
A. In the high frequency limit a; —>■ oo, the dominant contribution comes from 
stationary points of the phase, i.e. the points Xg for which 



-Xr\-\Xg-ij\- 



= 0 . 


The stationary points must then satisfy 


Xg — Xr 
\Xg - Xr 


Xg-y 

\xs - y| 


Thus by assumption |5.H there are no stationary points in the phase of the integral 


over the array A appearing in (23). Neglecting boundary effects, this integral goes 


to zero faster than any polynomial in w (see e.g. [2]). 
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We now show that in the high frequency limit w —>■ oo, we have ((xr,uj) —>■ 1/2. 
Recalling (211, we have 
(24) 

1 go{Xr,Uj)*p{Xr,Uj) - p{Xr,Uj)*go{Xr,Uj) 

C{Xr,Uj) = - H- 


go{Xr,uj)*go{Xr,uj) 


+ \z- Xr\ - \Xs - Xj.\)) 


— ^ dz dXgC{Xs)e-xp {iwCq'^ {\xs — Xr\ —\xs — z\ — \z — Xj.\)) , 
Cq J Ja 

where C{xg) collects geometric spreading terms and which is actually 

independent of the frequency uj. By assumption 5.1 the integrals over ^ in (24) do 
not have any stationary points. Thus if we neglect boundary terms, these integrals 
must go to zero faster than any polynomial in lo (see e.g. 0), meaning that 
((xr,uj) 1/2 as w —)■ oo. Thus TKMiCflo, w](y) —)■ 0 as w —)■ oo. □ 


6 . Autocorrelation measurements 

Up to this point we have assumed deterministic control over the source illumi¬ 
nations. In this section we relax this control by driving the array with stochastic 
signals. We start in section [6T] by recalling an ergodicity result of Gamier and Pa¬ 
panicolaou [To] which guarantees that if Gaussian stochastic processes are used to 
drive the sources, the realization average of the total field can be well approximated 
by time averages of the total field. Then in section [6^ we adapt the source pair illu¬ 
mination strategy to pairs of sources driven by two correlated Gaussian processes, 
with (known) correlation identical for different pairs. From these pairwise illumina¬ 
tions we measure empirical autocorrelations to obtain intensity measurements that 
are essentially (up to ergodic averaging) the same as those using the deterministic 
strategy of section [XT] 

6.1. Stochastic array illuminations. We consider array illuminations f(t) G 
given by a stationary Gaussian process with mean zero and with correlation the 
N X N matrix function 

(25) R{T) = (f{t)f{t + T)). 

Here (•) denotes the expectation with respect to realizations of /, and in an abuse 
of notation we have denoted by f{t) the time domain vector of signals driving 
the array. Since Rs.s'('r) = + =^(f+ T)fs(t)) = Rs/,s(-t) for 

s, s' = 1,..., A^, we have R(r) = R*(—r) and so R(a;) is a Hermitian N x N matrix. 

The total field u at the receiver arising from the array illumination f is, in the 
time domain, 

N 

(26) u{xr, t) = 

s=l 

where G is the Born approximation of the inhomogeneous Green function, i.e. 
G{xr,Xs,t) = ^ J doje 


GoiXr,Xs,Uj)+k‘^ / dzp{z)Go{Xr,Z,U!)Goiz,Xs,Uj) 


dt'G{Xr,Xs,t - 
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The empirical autocorrelation of u is 


(27) 


1 _ 

= — y u{Xr,t)u{Xr,t + T)dt, 


where T is a known measurement time. Following Gamier and Papanicolaou m, 
we formulate proposition |6.1| regarding the statistical stability and ergodicity of 


(271. This proposition is essentially the same as [ini Proposition 4.1], but we make 
small modifications to allow for complex fields and more general correlations in 
space. We include it here for the sake of completeness and the proof can be found 
in appendix [Xj 


Proposition 6.1. Assume f satisfies (251. The expectation (w.r.t. realizations of 
f) of the empirical autocorrelation (27) is independent of measurement time T: 

(28) 

where 


{ll}{Xr,r)) = 4'(®r,T), 


N 

(29) s,s'=i 

-n 


dt' 


duie 


dt"G{Xr, Xr, —t')G{Xr, Xs',T — t'')Rs^s'{t" — t') 
g{Xr,uj)*R{u:)g{Xr,u:). 


Furthermore, (27) is ergodic, i. e. 
(30) 


fj{Xr,T) 


T—>-oo 


> 4'(®r,T). 


6.2. Pairwise stochastic illuminations. We make Np illuminations each cor¬ 
responding to using only two distinct sources {i(m),jIrn)) G {l,...,iV}^, m = 
1,..., Np. The correlation matrix for the m—th experiment has the form 

(31) Rm.(w) = FmC(w)F)^, 

where F^ = [ei(Tn)T G and C{uj) is a known 2x2 Hermitian positive 

semidefinite matrix that represents the correlation between the two sources and is 
assumed to be the same for all experiments. For instance, if we send the same signal 
with power spectrum F(u}) from both sources in a pair, this correlation matrix is 


C(a;) = F{u!) 


By the ergodicity (301 of proposition |6.1[ when we measure the empirical au¬ 
tocorrelation il>rn of Ura at the receiver x^ for long enough time T, the empirical 
autocorrelation is close to an intensity measurement, i.e. 

(32) 'i>jn{Xr,Uj) = g{Xr,U})*FmC{uj)Fl,^g{Xr,U}). 

By using appropriate single source illuminations driven by a signal with known 
correlation, it is possible to measure 

(33) '^i{Xr,(^) = g*iXr,^)eieJ giXr,i^)- ioT i = I,..., N. 


From (32) and (33) we obtain the m—th measurement 


( 34 ) 


drn{Xr,0j) = ^m{Xr,0j) - Cn (w)$ (x^, w) - C 22 (w)§°(„) (x^ , w) 
= g{xr,uj)*F^D{oj)Fl,^g{xr,uj), 


















IMAGING WITH POWER GONTROLLED SOURCE PAIRS 


13 


where the matrix D is 2 x 2, Hermitian with zero diagonal, i.e. precisely of the same 


form as the matrix D we encountered in the intensity measurements case (11). 


Proceeding analogously as in section |3.1| and recalling that g = go + P we have 

dm{xr,uj) = {go + p)* FmD{uj)Fl^{go + p). 

Collecting the measurements for m = 1, • ■ •, Np and neglecting the quadratic term 
in p we have the approximate data 


(35) 


di{Xr,Uj) 

d2{Xr,Uj) 


{Xrj 


d{xr,oj) = M(i;r,a;) 


Refgo + 2p 
Im(go + 2p 


where the matrix M £ jg agajjj given by (131. Thus, the data (35) obtained 

by measuring the empirical autocorrelation (27) and using correlated pair illumina¬ 


tions, is essentially the same as the data obtained using deterministic source pairs 
(12). Hence the analysis of the matrix M of holds and we can use Kirchhoff 
migration as we did in for the intensity measurements case. 

Remark 6.2 (Uncorrelated background illumination). The proposed illumination 
strategy is robust with respect to noise and even allows to send the same Gaussian 
signal from the m—th source pair {i{m), j{m)) and independent Gaussian signals 
from all remaining sources on the array. If the independent signals have the same 
spectral density F{uj) as the source pair signal, the correlation matrix for the m—th 
experiment is 


(36) 


Rm (w) — F{uj) ( I -|- F. 


0 1 
1 0 


where I is N x N identity matrix. By subtracting from the autocorrelation for 
the m—th experiment, the autocorrelation for a reference illumination that sends 
independent Gaussian signals with correlation matrix F{u})\, it is possible to obtain 
m—th measurement (34) with 


□ (w) = F{uj) 


0 1 
1 0 


7. Additive noise 

Here we discuss the effects of additive instrumental noise in autocorrelated mea¬ 
surements of the total field. The total field at Xr resulting from illuminating with 
the m—th pair and tainted with additive noise is Um{xr,t) + f.{t)- We assume the 
noise ^ is a stationary Gaussian process with mean zero and spectral density 


(37) 


:(a;) = exp 


471 


Here Ic represents the correlation time of the noise (i.e. S(t) = {f,{t)f{t + t)) « 0 
for T 3> Ic) and ujq is the central angular frequency of the noise. If the noise f is 
independent of the signals used to drive the source pairs, it can be shown using the 
techniques of appendix [A] that 

1 pT 

— J dt{um{Xr, t) + fit)) {umiXr, t + t) + f{t + t)) t) + S(t), 
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where is given by (291. 

Assuming the same form of instrumental noise in the single source reference 
measurements, the m—th measurement dmixr,uj) is 

dm{xr,u}) = (go +p)*FmD(a;)F^(go +p) +CE{oj). 

for some C € M. Neglecting the terms which are quadratic in p and going back to 
the time domain we have 






duje 


+9o(£r,w)*F^D(w)F((,p(£r,w) 

+p{Xr,Uj)*FrnD{uj)F'l^go{Xr,Uj) + CE{t), 


with the slight abuse of notation of using dm for both time and frequency domain 
quantities. The second and third terms in the integrand are incident-scattered held 
correlations and contain the available information about the scattering potential 

piy)- 

For simplicity, we now focus on the case where the source pair signals have 
correlation matrix 


D(u;) = F{ijj) 




Such correlation corresponds to sending a signal from one of the sources in a pair and 
a copy of the same signal delayed by (j) from the other source. For a point scatterer 
at y, the incident-scattered terms have peaks at delay times r(y) corresponding 
to differences between travel times of a rehected path and direct path, i.e. for the 
m—th experiment the peaks occur at the four possible delays 


r{y) 


±((l®l(m) -y\ + \y-Xr\ - \x,(^m) “ X r-\) / Cq + (j)), 
±((l®^(m) -y\ + \y-Xr\ - \Xj(^m) “ ®r|)/co “ </>). 


(38) 


Consider then the minimal delay time Tmin(y) given by 

Xs-y\ + \y- Xr\ - \x^, - Xr 


i(y) = ^ min 

Xs,X^i£A 


’ zb ( 


Co 


that is the minimal delay time we expect the incident-scattered correlations to peak. 
If we assume the additive noise decorrelates much faster than the first incident- 
scattered arrival from y (i.e. Ic ^ 'imin(y)), then the information of the scatterer 
p(y) contained in dmixr,T) is essentially unchanged (up to ergodic averaging). 
Hence we can stably image using the proposed method at y provided Tmin(y) ^ 4- 


8 . Numerical experiments 

Here we include 2D numerical experiments of our proposed imaging routine for 
scalings corresponding to acoustics ()8.1) and optics ()8.2). We demonstrate the 
stochastic source pair illumination strategy for the acoustic regime, i.e. we compute 
the autocorrelations for time domain data. In the optic regime this is an expensive 
calculation, so we use instead power spectra (i.e. deterministic illuminations). 
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8.1. Acoustic regime. For imaging in an acoustic regime, our choice of physical 
parameters corresponds to ultrasound in water. We choose the background wave 
velocity to be cq = 1500 m/s. The central frequency for all signals (sources and 
additive noise) is 3 MHz, which gives a central wavelength of Aq = 0.5 mm. We 
center a source array A at the origin consisting of 41 sources at coordinates Xg = 
(0, —lOAo + (s — l)Ao/2) for s = 1,..., 41. A single receiver is located at the 
coordinate x^ = (— 20Ao, — 20Ao) (see figure]^. 

We generate a stationary Gaussian time signal f(t) with mean zero and correla¬ 
tion function 

F(t) = exp - 

using the Wiener-Khinchin theorem. The correlation time tc ~ 1.25 which gives 
the signal an effective frequency band [1, 5] MHz. We generate time signals of length 
2T for T « 260 /iS with 8001 uniformly spaced samples. This sampling is enough 
to resolve the frequencies in the angular frequency band B, while T is sufficient to 
observe ergodic averaging (see (j^. By placing the same realization of this signal 
f{uj) at the locations iqm) and alj(m) "we generate the pair illumination fm (oj) = 
f(oj)(ei(m) + Similarly, by placing an independent realization of /(w) at 

location Xi we generate the single source reference illumination /f(w) = f{Lo)ei. 

For all experiments, synthetic data is generated in the frequency domain using 
the Born approximation. We assume 3D wave propagation for simplicity so that Go 
is given by ([^ for d = 3. The m—th measurement is obtained through the formula 

d^{xr,u:) = ^m{xr,uj) - $°(^)(xr,a;) - $°(„)(®^,w), 

where 



-p ^ 

{gQ{Xr,Uj)+p{Xr,aj)) fm{uj) 


T ^ 

{go{Xr,Uj) +p{Xr,Uj)) f°{u}) 


with Qq and p defined by (H and (§ respectively. 

For these simulations we use the full set of pair illuminations, which for A = 41 
source locations, generates a measurement matrix M(;Er, w) G ]R® 20 x 82 ^ 

Moore Penrose pseudoinverse to recover p + Cgo for each uj G B. When the num¬ 
ber of sources N and thus the dimension of M is large (recall M S ^N{N-1)/2XN^^ 
the pseudoinverse could become computationally expensive. However, the system 
is sparse as it contains only 4 non-zero elements per row, so linear least square 
solvers that exploit sparsity (e.g. CGLS m) may be more efficient than our ap¬ 
proach. Furthermore, as discussed in (j^we can reduce the size of M to 2N x 2N 
while keeping the nullspace of M one-dimensional by using an appropriate subset 
of source pairs. 

We form an image aty GW = {(100Ao-l-iAo/2.5, jAo/2.5), for i,j = —25,..., 25} 
using the Kirchhoff migration functional (( |2.2[ ), summed over the bandwidth band 

B, 

rKM[p + C9o](y) = / dwFKM[p + C9 o,w](y). 

JB 

For our first experiment, we place a single point reflector at the location y = 
(100Ao,0), with refractive index perturbation p{y) = 1 x 10“®. The migrated 
image (figure [5}r) indeed exhibits the cross-range (Rayleigh) resolution estimate 
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XoL/a « 5Ao and range resolution estimate co/|/B| « IAq. Note that there is a 
trade-off in the choice of the reflectivity: p has to be sufficiently small so that the 
quadratic terms in p can be neglected in (12). However the smaller p is, the longer 
the acquisition time T has to be in order to better observe the refiected-incident 
correlations in the data. 

In our second experiment (figure]^), we consider two oblique reflectors located 
at Vi = (99Ao,—2Ao) and 1/2 = (103Ao,4Ao) each with p{yi) = 1 x 10“®. We 
include a reconstruction of an extended scatterer (line segment) in figure Here 
the line segment is generated as a set of point reflectors each with p{yi) = 1 x 10“® 
uniformly spaced by Ao/8. 


rKM[p](y) rKM[p + CsoKy) 



90 95 100 105 110 90 95 100 105 110 



90 95 100 105 110 90 95 100 105 110 


Figure 5. Kirchhoff images of (a) one point and (b) two point re¬ 
flectors, whose true positions are indicated with crosses. The left 
column uses the full waveform data p, while the right column use 
the recovered data p + (qq. The horizontal and vertical axes dis¬ 
play the range and cross-range respectively, with scales in central 
wavelengths Aq. 


We now demonstrate the robustness of the proposed method with respect to 
additive noise (see section . Here we have taken a realization of the data for 
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rKM[p](27) rKM[p +Cfifo](y) 



90 95 100 105 110 90 95 100 105 110 


Figure 6. Kirchhoff images of an extended reflector. The left 
column uses the full waveform data p, while the right column use 
the recovered data p + Cgo- The horizontal and vertical axes dis¬ 
play the range and cross-range respectively, with scales in central 
wavelengths Aq. 


a single point reflector (c.f. figure [^) and perturbed each measurement with 
additive noise as follows. The m—th signal Umixr,to) has total power pm = 
f \um{xr,ui)\‘^duj. We construct a Gaussian signal ■fm(i) with mean zero, spectral 
density (371, 4 ~ 1-25 ps and total power 1. This allows to obtain the perturbed 
total field Um(xr,^) + for some V > 0. The to— th measurement with 

additive noise is thus dm{xr,^) = |wm(®r,w)p + ■ Thus the ratio of the 

signal power to the noise power is 1/u. The signal-to-noise ratio (SNR) is then 


SNR,„ = -101ogio(u)dB. 


Figure [7] shows the reconstruction from data with SNR™ = 0 dB for each to, 
meaning that the signal and the noise have the same power. 

Lastly we perform an experiment that sends as the to— th illumination the usual 
correlated pair illumination fm, and uncorrelated noise from the remaining sources 
on the array A (see remark 6.2). To generate this illumination we place the same 
realization of the signal /(w) at the locations and Xj(^^p and independent 


realizations of /(w) at the remaining source locations. Similarly, a reference illumi¬ 
nation is generated by placing independent realizations of /(w) at all locations on 
the array A. By measuring the autocorrelation of the resulting fields we obtain data 
that is essentially the same form as dm{x^, w). Figure|^shows this experiment with 
the single point reflector located ai y = (IOOAq, 0) and reflectivity p{y) = 1 x 10“®. 


8.2. Optic regime. For imaging in an optic regime, we use the background wave 
velocity cq = 3 x 10® m/s and central frequency « 589 THz which gives a central 
wavelength Aq ~ 509 nm. Our source array A is again centered at the origin, but 
now consists of 1001 sources located at coordinates x^ = (0, —500Ao + (s — l)Ao) 
for s = 1,..., 1001, and we set = (—IOOOAq, —IOOOAq). 
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rKM[p](27) rKM[p +Cfifo](y) 



90 95 100 105 110 90 95 100 105 110 


Figure 7. Additive noise: (left) array response vector migra¬ 
tion rKM[p](y), (right) recovered array response vector migration 
rKM[p + Cgo]{y) for SNR^ = OdB. The horizontal and vertical 
axes display the range and cross-range respectively measured in 
central wavelengths Aq. 


rKM[p](y) ^km[p +Cgo]{y) 



90 95 100 105 110 90 95 100 105 110 


Figure 8. Uncorrelated background illumination: (left) array re¬ 
sponse vector migration rKM[p](y)) (b) recovered array response 
vector migration rKM[p + Cgo]{y) for SNR^ = OdB. The horizon¬ 
tal and vertical axes display the range and cross-range respectively 
measured in central wavelengths Aq. 


We generate intensity data d{xr,uj) as 


— (fl^O T p) T ^t(m)) (PO T p) ^i(m) (PO T p) 


°j (m) 


for 100 (angular) frequencies oj uniformly spaced in the frequency band [429, 750] 
THz. This corresponds to performing the source pair experiments (source pair 
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illuminations and single source reference illuminations) for 100 different monochro¬ 
matic visible light sources with wavelengths A € [400, 700] nm, equally spaced in 
frequency. Since there are a large number of sources in this setup (N = 1001), 
we implement the strategy discussed in Q to reduce the number of source pair 
experiments from Np = N{N — l)/2 to Np = 2N. 

As before, we use the pseudoinverse to recover p + Cgo for each frequency oj € 
B, and then use the Kirchhoff migration functional ([ |2.2[ ) to form an image. Here 
we use the image window W = {(5000Ao-I-jAo/2.5, jAo/2.5), for i,j = —25,..., 25}. 
In figure l^b) we demonstrate the migrated image for two point reflectors placed 
at Vi = (4098Ao,3Ao) and y 2 = (5004Ao, —5Ao) each with reflectivity p{yi) = 1 x 
10“^^. Although we are significantly undersampling the data in frequency and the 
source spacing is larger than Ao/2, the spot sizes still exhibit the Kirchhoff migration 
resolution estimates ([ 2.2 ) of XoL/a « 5Ao in cross-range and co/|H| « 2Ao in range. 


10 


-5 


-10 


rKM[p](y) 


^km[p + Cgo]{y) 



Figure 9. Optic regime: (left) array response vector migration 
rKM[p](y)) (b) recovered array response vector migration rKM[p + 
Cgoliv)- The horizontal and vertical axes display the range and 
cross-range respectively measured in central wavelengths Aq. 


9. Discussion 


By sending correlated signals from different pairs of locations we have shown 
that from intensity data we can recover full waveform data by solving a linear 
system. This linear system has a known one-dimensional nullspace provided the 


sources and receiver satisfy the distance conditions given by assumption 4.1 which 
allows for the recovery oip +(go- We show this quantity is enough to use standard 
migration techniques (e.g. Kirchhoff migration Tkm) provided the sources and 
receiver satisfy the additional geometric conditions of assumption |5.1| Thus we 
obtain full waveform resolution estimates for an image formed from intensity-only 
data. 

Our method relies only on knowledge of paired source locations and the correla¬ 
tion of the signals being sent. This allows us to relax illumination control by using 
paired stochastic signals. By measuring autocorrelations of the resulting fields, we 
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obtain essentially the same intensity data as with using deterministic source pairs. 
These stochastic illuminations can be created e.g. by using a configurable mask 
that is parallel to the wave fronts of an incoherent plane wave. 

The linear system we solve has size 2N x 2N and is very sparse (up to 4 non¬ 
zero entries per row). In our simulations we used M^, however sparse solvers such 
as CGLS (see e.g. [HI) could be used. To form the system we need at least 
3iV different illuminations, 2N pair illuminations plus N reference illuminations. 
However, in our illumination strategy, the phase of the source signals does not need 
to be known. We replace the direct phase control by the natural phase modulation 
that comes from the different positions of the signals. 

We use the geometric imaging conditions (assumption |5.1[ ) to show the nullspace 
of M does not affect imaging via Tkm- This assumption imposes some restrictions 
on the juxtaposition of the sources and receiver and in turn on the forms of illumi¬ 
nations we can consider. For example, using a stationary phase argument, it can be 
shown the autocorrelation of the total field is negligible if spatially continuous array 
illuminations (rather than paired point sources) are used. In future work, we would 
like to address this more thoroughly to determine if more general illuminations can 
be used. It may also be interesting to see if the source pair strategy we propose 
will work for other imaging setups. 
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Appendix A. Proof of proposition 16.11 

In this appendix we prove proposition |6.1| which details the statistical stability of 
the measured autocorrelation (27) with respect to realizations of the illumination 
/. The theorem and proof are patterned after the result by Gamier and Papanico¬ 
laou [ini Proposition 4.1], only we make small modifications to allow for complex 
fields and the form (25) of the correlation function R(t). 


Proof. Since we are assuming / is a stationary process in t, the resulting total field 
u is also a stationary random process in t. So we have 

{u{Xr, t)u{Xr, t + t)) = {u{Xr, 0)u{Xr, t )), 

which allows us to compute 


1 F 

{lp{Xr,T)) = ^ J ^dt{u{Xr,t)u{Xr,t + T)) 

1 F _ 

= ^y dt{u{Xr,0)u{Xr,T)) = {u{Xr,0)u{Xr,T)). 
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So (27) is independent of T. By expressing the quantity {u{xr, 0)u{xr, t)) through 
the Green’s function G we verify (28): 


{tp{Xr,T)) = ^ J dt' J dt''GiXr,Xp,-t')G{Xr,Xp,,T 

P,P' — 1 ^ 

= ^ J dt' J dt”G{Xr,Xp,—t')G{Xr,Xpi,T — t”)Rpy{t” — t') 

p,p' — l 


= ^ / due 


To show the ergodicity (30), we need to compute the variance of ijj. We first 
compute the covariance as 
(39) 

Cov('!/:(a?r, r), V'(®r, T + At)) = ^ / / dtdt' / dsds'dudu' 

p^p' ^q^q—\ 

X G{Xr, Xp, s)G{Xr, Xpi,U — T)G{Xr, Xq, s')G{Xr, Xgf ,u' — T — At) 

X - S)jp>{t - U)fqit' - S')fq.{t' - U')) 

- {fp{t - S)7p'(i - U))(jq{t' - S)fq>{t' - U 

The product of the second order moments is 

{fpit - S)]p,it - U)){fq{t' - S')fq>{t' - U')) = Rp'.piu - s)Rq^q'{s' - u'). 

Since f{t) is Gaussian (in time), the fourth order moment is given by the complex 
Gaussian moment theorem (see e.g. [22]) as 

{fp{i - s)7p>{t - - s')fq^{t' - u')) = Rp>^p{u - s)Rq^q>{s' - u') 

“h Rq p(t — t — S “h S ^Rp'^q'(t — t — U ti). 


We now integrate over the t, t' variables to obtain 


/^^^'((/p7-'S)/p'7-w)/9(i'-s')/g'(^'-w')) 

- {fp{t - s)fp>it - u))ljq{t' - s')fq,(t' - u' 


= J dt j dt'Rq^p{t - t' - S + s')Rp'^q> {t' -t - u' + u) 

= j dio J duj' sinc^iuj - uj')T)e^‘^'^^-^''>e-^‘^^^-^'^Rq,p{uj)Rp>y{uj'). 






22 


PATRICK BARDSLEY AND FERNANDO GUEVARA VASQUEZ 


Plugging this into (39) we obtain 


Cov {^ijj{xr,T),tp{xr,T + At)'^ = ^ JdujJ dw'sinc^ ((w — w')T) 

p,p',q:q'=i 


X G{Xr, Xp, U!')G{Xr, Xpi , U})G{Xr, Xq,U:)G{Xr, Xqi , Uj)Rq^p{uj)Rpi ^q' 

= J duj J dcj'sinc^ ((w — w')T) (^g{xr,uj)*R{oj)g{xr,oj')) 

X [g{xr,u}yR{uj')g{xr,u}'))G‘^^'" , 

where g = go + p is given by ([^ and (|^, and (R(a;)).^. = Rij{uj) is a 
Hermitian matrix for each to. Then taking T —>■ oo we compute the variance as 

r Var t)) ——^ / dw g(£r,a;)*R(a;)g(£r, w) 


and so the variance is 0(1/T) as T —>-00. This establishes (30). 


□ 
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